In [2]: 
import scipy.signal
In [8]: 
def plotspec(x, Ts):
    fig = figure()
    ax1 = fig.add_subplot(211)
    ax1.plot(x)
    
    q = fft.fft(x)
    ax2 = fig.add_subplot(212)
    ax2.plot(fft.fftfreq(len(x), Ts), abs(q))
    
time = 3.0
Ts=1.0/10000.0
x = random.uniform(-1.0,1.0, time/Ts)
plotspec(x, Ts)

figure()
b = scipy.signal.remez(100, [0,0.1,0.105,0.5],[1,0])
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)


figure()
b = scipy.signal.remez(100, [0,0.12,0.13,0.25,0.255, 0.5],[0,1,0])
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)


figure()
b = scipy.signal.remez(100, [0,0.37,0.38,0.5],[0,1])
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)

Notes that remez() takes slightly different arguments in scipy.signal than it does in MATLAB:

  • Top end of frequency specification is 0.5 instead of 1.0 (ie: Nyquist frequency), so all those terms are halved
  • The third argument takes half as many inputs as the previous one: it expects that we're defining levels between the bands, or something

3.8. Mimic the code in filternoise to create a filter that

  1. passes all frequencies above 500Hz
  2. passes all frequencies below 3000Hz
  3. rejects all frequences between 1500 and 2500 Hz
In [20]: 
def plotfilter(b, Ts):
    figure()
    time = 3.0
    x = random.uniform(-1.0,1.0, time/Ts)
    y = scipy.signal.lfilter(b,1,x)
    plotspec(y, Ts)


def genlpf(ntaps, topf, Ts):
    q = topf/(1/Ts)
    b = scipy.signal.remez(ntaps, [0, q, 1.08*q, 0.5], [1,0])
    return b

def genhpf(ntaps, bottomf, Ts):
    q = bottomf/(1/Ts)
    b = scipy.signal.remez(ntaps, [0, q, 1.08*q, 0.5], [0,1])
    return b

def genbandpassf(ntaps, bottomf, topf, Ts):
    q = bottomf/(1/Ts)
    r = topf/(1/Ts)

    b = scipy.signal.remez(ntaps, [0,q*0.98,q,r,r*1.02,0.5], [0,1,0])
    return b
In [23]: 
Ts = 1.0/10000.0

plotfilter( genhpf(100,  500.0, Ts), Ts)
plotfilter( genlpf(100, 3000.0, Ts), Ts)
plotfilter( genbandpassf(100, 1500, 2500, Ts), Ts)

3.9. As 3.8, but with Ts=1/20000.0

In [26]: 
Ts = 1.0/20000.0

plotfilter( genhpf(100,  500.0, Ts), Ts)
plotfilter( genlpf(100, 3000.0, Ts), Ts)
plotfilter( genbandpassf(100, 1500, 2500, Ts), Ts)

3.10. Let x1 be a cos of f=800, x2(t) of f=2000 and x3 of f=4500. Let x(t) = x1 + 0.5x2 + 2x3

In [27]: 
def plotfilter2(b, Ts):
    figure()
    time = 3.0
    tp = 2*pi*linspace(0, time, time/Ts)
    x = cos(800*tp) + 0.5*cos(2000*tp) + 2*cos(4500*tp)
    y = scipy.signal.lfilter(b,1,x)
    plotspec(y, Ts)
    
plotfilter2( genhpf(100,  500.0, Ts), Ts)
plotfilter2( genlpf(100, 3000.0, Ts), Ts)
plotfilter2( genbandpassf(100, 1500, 2500, Ts), Ts)